#include <iostream>
#include <cuda.h>

__global__ void MatrixAdd_d(float *A, float *B, float *C, int N)
{
	int i = blockIdx.x*blockDim.x + threadIdx.x;
	int j = blockIdx.y*blockDim.y + threadIdx.y;
	int index = i*N + j;
	//printf("N=%d i=%d j=%d index=%d\n",N,i,j,index);
	if(i<N && j<N)
	{
		C[index] = A[index] + B[index];
		//printf("A[%d]=%f B[%d]=%f C[%d]=%f\n",index,A[index],index,B[index],index,C[index]);
	}
}

__global__ void smoother(float *A, float *B, int N, int kernel_width)
{
	int row = blockIdx.x*blockDim.x + threadIdx.x;
	int col = blockIdx.y*blockDim.y + threadIdx.y;
	int index = row*N + col;
	unsigned int window_width = (2*kernel_width+1);
	printf("N=%d i=%d j=%d index=%d width=%d\n",N,col,row,index,window_width);
	if(row<N && col<N)
	{
		float value = 0;
		int num_kernel_cells = window_width * window_width;
		for ( unsigned int i=0; i<window_width; i++ )
		{
			int kernel_row = (row + i) * N;
			for ( unsigned int j=0; j<window_width; j++ )
			{
				int kernel_col = (col + j);
				value += A[kernel_row + kernel_col];
			}
		}
		B[index] = value / ((float)num_kernel_cells);
		printf("value=%f\n",B[index]);
	}
}

/*
 * Function: InitField
 *
 *  Initiliaze the input field to well known values for debugging
 *  and visualization purposes.
 */
__host__ void InitField ( unsigned dim, float* pField )
{
  for ( unsigned i=0; i<dim; i++ )
  {
    for ( unsigned j=0; j<dim; j++ )
    {
      pField[i*dim+j] = j + i;
    }
  }
}

       
int main()
{
	float *old_h, *new_h; // pointers to host memory; a.k.a. CPU
	float *old_d, *new_d; // pointers to device memory; a.k.a. GPU
	int blocksize=16, kernel_width = 8, n=2064, smoothed_n = n - (2*kernel_width), i, j, index;
	
	// allocate arrays on host
	old_h = (float *)malloc(sizeof(float)*n*n);
	new_h = (float *)malloc(sizeof(float)*smoothed_n*smoothed_n);
	
	// allocate arrays on device
	cudaMalloc((void **)&old_d,n*n*sizeof(float));
	cudaMalloc((void **)&new_d,smoothed_n*smoothed_n*sizeof(float));
	dim3 dimBlock( blocksize, blocksize );
	dim3 dimGrid( smoothed_n/dimBlock.x, smoothed_n/dimBlock.y );
	
	InitField(n, old_h);
	InitField(smoothed_n, new_h);

	// copy and run the code on the device
	cudaMemcpy(old_d,old_h,n*n*sizeof(float),cudaMemcpyHostToDevice);
	cudaMemcpy(new_d,new_h,smoothed_n*smoothed_n*sizeof(float),cudaMemcpyHostToDevice);
	smoother<<<dimGrid, dimBlock>>>(old_d,new_d,smoothed_n,kernel_width);
	cudaThreadSynchronize();
	cudaMemcpy(new_h,new_d,smoothed_n*smoothed_n*sizeof(float),cudaMemcpyDeviceToHost);

	// print out the answer
	for(i=0;i<smoothed_n;i++)
	{
		for(j=0;j<smoothed_n;j++)
		{
			index = i*n+j;
			printf("old_h[%d]=%f new_h[%d]=%f\n",index,old_h[index],index,new_h[index]);
		}
	}
	
	// cleanup...
	free(old_h);
	free(new_h);
	cudaFree(old_d);
	cudaFree(new_d);
	return(0);
}